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Dynamic turbulence modelling in large-eddy 
simulations of the cloud-topped atmospheric 

boundary layer 

By M. P. Kirkpatrick, N. N. Mansour, A. S. Ackerman f, and D. E. Stevens $ 


1. Motivation and Objectives 

The use of large eddy simulation, or LES, to study the atmospheric boundary layer 
dates back to the early 1970s when Deardorff (1972) used a three-dimensional simulation 
to determine velocity and temperature scales in the convective boundary layer. In 1974 
he applied LES to the problem of mixing layer entrainment (Deardorff 1974) and in 1980 
to the cloud-topped boundary layer (Deardorff 1980 b). Since that time the LES approach 
has been applied to atmospheric boundary layer problems by numerous authors. 

While LES has been shown to be relatively robust for simple cases such as a clear, 
convective boundary layer (Mason 1989), simulation of the cloud-topped boundary layer 
has proved more of a challenge. The combination of small length scales and anisotropic 
turbulence coupled with cloud microphysics and radiation effects places a heavy bur- 
den on the turbulence model, especially in the cloud-top region. Consequently, over the 
past few decades considerable effort has been devoted to developing turbulence models 
that are better able to parameterise these processes. Much of this work has involved 
taking parameterisations developed for neutral boundary layers and deriving corrections 
to account for buoyancy effects associated with the background stratification and local 
buoyancy sources due to radiative and latent heat transfer within the cloud (see Lilly 
1962; Deardorff 1980a; Mason 1989; MacVean & Mason 1990, for example). In this pa- 
per we hope to contribute to this effort by presenting a number of turbulence models 
in which the model coefficients are calculated dynamically during the simulation rather 
than being prescribed a priori. 

The dynamic procedure, first proposed by Germano et al. (1991), computes model 
coefficients using information contained in the resolved velocity and scalar fields. In this 
sense dynamic models can be considered self-calibrating, a feature that makes them an 
appealing choice for dealing with the complex interactions between the hydrodynamics, 
radiation and cloud microphysics occurring within clouds. Nevertheless, while dynamic 
turbulence models have been used with considerable success for complex engineering 
flows (see Boivin et al. 2000; Branley & Jones 2001, for example), their application to 
atmospheric flows has been very limited. 

The objective of this paper is to present a set of dynamic turbulence models written 
in a form appropriate for large eddy simulations of the atmospheric boundary layer in 
which the flow is described using equations of motion in anelastic form . The models are 
tested using simulations of a nocturnal marine stratocumulus cloud observed during dur- 
ing the second Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) field 
experiment. This test case includes a number of features that typically cause problems 
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in simulations of the atmospheric boundary layer, including the presence of a strong 
temperature inversion at the level of the cloud top, and positive feedback loops involv- 
ing turbulent entrainment across the inversion, radiation and cloud microphysics. As 
such, these simulations provide a good test for the utility of our models. Simulations 
are performed using the atmospheric boundary layer LES code, DHARMA, (Stevens & 
Bretherton 1999; Stevens et al. 2000), with dynamic turbulence models computed using 
the new LLAMA code presented here. 


2. Dynamic turbulence models for the anelastic equations 

The dynamics of the cloud-topped atmospheric boundary layer can be described using 
equations for conservation of mass, momentum, liquid water potential temperature and 
total water mixing ratio. These equations are written in the anelastic form of Ogura & 
Phillips (1962) in which the thermodynamic variables such as pressure p are decomposed 
into an isentropic base state pp (corresponding to a uniform potential temperature 9p) 
and a dynamic component. Following Clark (1979), the dynamic component is further 
decomposed into an initial environmental deviation in hydrostatic balance p\ and a time- 
evolving dynamic perturbation p 2 to give 

p(x,y,z,t ) = po(z) + pi{z) + p 2 (x,y, z,t). ( 2 . 1 ) 


After subtracting hydrostatic balances, the resulting continuous equations are written, 


dg 0 Uj djeoUjUj) 

dt dxj 

dgpOi d(QpO*Uj) 

dt dxj 

dgpqt d(g 0 q t Uj) 

dt dxj 

d(gpuj) 

dxj 


dp2 

dxi 


+ 


gp0v2 


= Her 


= H„ 


= 0. 


tijkqp f j^k T 


( 2 . 2 ) 

(2.3) 

(2.4) 

(2.5) 


Here w, are the Cartesian components of the velocity vector, g the density, Sij the Kro- 
necker delta, the permutation tensor, fj the Coriolis parameter, g the acceleration 
due to gravity, qt the total water mixing ratio and 0 ; * = ( 6i — dp) /6p a scaled liquid wa- 
ter potential temperature. Total water mixing ratio is the sum of the liquid and vapour 
mixing ratios, 


g c T 9v 

qt = q c + q v = 


( 2 . 6 ) 


Qd 

where g c , g v and qd are the density of the condensed water, water vapour and dry air 
respectively. Liquid water potential temperature is defined as 


0i 


Lq c 

CpdTTp 


(2.7) 


Here L is the latent heat of vaporisation, C p d the specific heat at constant pressure for 

dry air and np = ( jp 8 —) P with p re f a reference pressure and R,i the gas constant of 
dry air. The virtual potential temperature 0 V appearing in the buoyancy term of the 
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momentum equations is given by 


0 V — 9 + 9q 



( 2 . 8 ) 


where R v is the gas constant for water vapour. The terms H Ui , Hg * and H qt are source 
terms that include parameterisations for physical processes such as radiation, precipita- 
tion and the effects of the global atmospheric circulation. 

In a large eddy simulation, the equations are filtered to remove from the solution those 
fluctuations that cannot be resolved by the numerical method. For the anelastic equations 
it is convenient to use a density- weighted or Favre filter, where a Favre filtered variable 
is defined as <f> = Q<t>/6- Application of this filter to the equations gives 


dg 0 Ui d(g 0 UiUj) 


dt 

dt, 

dg 0 qt 
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dp 2 , x Ihflvi 
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d(QpUj) 
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with subgrid scale stresses and fluxes given by 


Tij £>q ( Ill'll j , 

= Qo {°i u j ~ tfuj) » 

7 q t = Qo {otUj - QtUj ) . 


(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

(2.12) 


(2.13) 

(2.14) 

(2.15) 


We have assumed here that fj is constant so that the Coriolis term is linear - an assump- 
tion that is valid when the domain is small compared to the Rossby radius of deformation. 
The source terms H are shown without an overbar since they are parameterisations rather 
than exact terms. 

The Smagorinsky model (Smagorinsky 1963) for the subgrid scale stress is 


X j,j — Tij d/jTgg — 2/0,, ErnJdt jCd a , 

where Dij is the strain rate tensor, 

n _ I ( i _ lx 
13 2 \ dxj dxi ) 3 13 d'Xg ' 


(2.16) 


(2.17) 


Here and elsewhere in this paper the superscript ° is used to denote the anisotropic part 
of a tensor. The eddy viscosity K rn is given by 


I< m = C A 2 


D 


(2.18) 


with 


D 


2DijDij and C a dimensionless coefficient. Stratification effects are ac- 
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counted for using the correction of Lilly (1962), 

C B = ( 1 - Riegs) 1 ' 2 Risge < 1 (2.19) 

C B = o Risgs > 1, (2.20) 


where the subgrid scale flux Richardson number is estimated using 


_ g/Oo dO v /dz 

Prsgs\D\2 ' 


( 2 . 21 ) 


Pr sgs is the subgrid scale turbulent Prandtl number, which we set to a constant value of 
Pr sgs = 0.4. 

The mixed model of Bardina et al. (1983), combines the Smagorinsky model with the 
scale similarity model of Bardina et al. (1980). The model is derived by rewriting the 
subgrid scale stress tensor in terms of resolved and unresolved parts of the velocity field 
Ui = Uj + u\ using the decomposition of Germano (1986). For the anelastic equations this 
decomposition takes the form 


r u - • <7*. • R[ r (2.22) 

where the modified Leonard, cross and SGS Reynolds terms are 

Rij — Qp /-// u.j ) (2.23) 

C*j = g 0 (uiu'j + u' i u j - Uiu'j - u'iUj) (2.24) 

R* ij =eoW^' j -u' i u' j ). (2.25) 


The modified Leonard term, which contains only resolved scale quantities, can be calcu- 
lated explicitly leaving the modified cross and SGS Reynolds terms to be parameterised 
using the Smagorinsky model, 


= L*a - 2g 0 CA 2 


D 


DijCs- 


(2.26) 


Martin et al. (2000) note that, with the mixed model, the isotropic part of the subgrid 
scale stress r,y can be represented at least partially by the isotropic part of the modified 
Leonard term. In fact they obtained the closest agreement with experiment for their 
compressible flow test cases when they assumed that the isotropic part of L*j models the 
whole of \5ijTkk- We have adopted the same approach in Eq.(2.26). 

In dynamic versions of these models, the model coefficient C is calculated dynamically. 
To achieve this, a test filter is applied to the velocity and scalar fields to extract informa- 
tion from the smallest resolved scales. Application of a spatial test filter (denoted here 
by a caret) to the filtered momentum equations gives 


dgyUj d{g 0 UjU 3 ) 

dt dxj 


dp* 2 

dxi 


+ Sag 


8p®v2 

do 


e ijkQpU k f] T Hu, 


d? jj 
dxj ' 


(2.27) 


In order to rewrite this equation in a form similar to Eq.(2.9) we adopt a Favre test filter 
(f>= g<t>/Q giving 


dg 0 Uj d(g 0 UjUj) 
dt dxj 


dpi 

dxi 


+ Si3g 


Qq@v2 


Ci j k g {) Uk f j “1 Hy , 


d%j_ 

dxj ' 


(2.28) 


Since density varies only in the 2 direction, this equation can be simplified considerably 
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by using two-dimensional horizontal test filter rather than a three-dimensional one. The 
Favre test filter is then equivalent to a spatial test filter, 


and the density field is unchanged by the test filtering operation, that is, 

Qo = Qo- 

From Eqs.(2.29) and (2.30), with horizontal test filtering Eq.(2.28) becomes 


dg 0 Ui d(g 0 UiUj) 


dt 


dxj 


Jj _ &V2 , r — — ~ f tt '-"ij 

— + °i3t>29 ~ e ijkQo u kfj + H Ui — . 


dfi , 


Extracting the resolved stress = g 0 (v,iUj — UiUj ) then gives 


dg 0 Uj d(g 0 UjUj) _ dp. 

dt dxj dxi 


dx.j 


dr-H dLi 


2 + S i3 g 2 g - c,. lk o it h k f, + 


(2.29) 

(2.30) 

(2.31) 

(2.32) 


Alternatively, application of the grid and test filters together to the continuous equations 
gives 


dg 0 Ui d(g 0 iiiUj) dp: 


'j <r 
dt 


dxj 


— 9 — e ijkQo u kfj + H Ui — 


dxj 


where 


h'ij 9(] ^UiUj UiUj^j . 


Comparison of Eqs.(2.32) and (2.33) gives the Germano identity, 


L ij — d 'i j Tj j . 


(2.33) 

(2.34) 

(2.35) 


The dynamic procedure assumes that the same parameterisation that is used to pa- 
rameterise the SGS stress at the grid filter level ry can be used to model the test level 
stress Tij. The test level stress for the Smagorinsky and mixed models are written as 


= ~2g 0 CA 2 


D 


DijC B 5 


and 


Tij - L*J - 2g 0 CA 

respectively. Here Cb is given by 

C B =(l- Ri 
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DijC B 5 


‘‘sgs 


1/2 
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Cb = 0 Ri S gs > 1 ; 

with the test level Richardson number given by 

9 /do dd v /dz 


Fti — 

ll ' t sgs — 


Pr. g ,\D \ 2 


(2.36) 

(2.37) 

(2.38) 

(2.39) 

(2.40) 


The test level modified Leonard term in the mixed model L*J is found by decomposing 
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in a manner similar to the decomposition of Tij (see Eq.(2.22)). Substituting u t = Ui + u" 
into Eq.(2.34) gives 

Tn = L*J + C*J + R*J, (2.41) 

where the modified Leonard, cross and Reynolds terms at the test level are 


L'ij Qo(uiUj UiUj) (2.42) 

C\J = 0o(uiu'j + u"uj — Ui.u " — u"uj) (2.43) 

= (2-44) 

In Eqs. (2.36-2.44) we have again used the fact that we assume a horizontal test filter so 
that Eqs. (2.29-2.30) can be used to simplify the equations. 

The decomposition shown in Eqs. (2.42-2.44) is different from that derived by Zang 
et al. (1993), who used a decomposition based on rq = tq + v! i at both the grid and test 
level, rather than the form zq = zq + u" used here for the test level. As noted by Vreman 
et al. (1994), the latter is more consistent with the dynamic procedure since the test level 
stress should be written in terms of test filtered velocity components in order to make 
the test level parameterisations analogous to those used at the grid level. 

Substituting the parameterisations for the subgrid scale stress Tq (Eq.(2.26)) and test 
level stress Tq (Eq.(2.37)) into the Germano identity (Eq.(2.35)) gives 


Qo ( A? 



Qo(uiUj - UiUj) - 2g 0 CA 2 
-Q 0 (uiUj - UiUj) + 2g 0 CA 2 


d,,,c b 

DijCs- 


(2.45) 


Removing the constant factor g 0 , and contracting tensors using the least squares approach 
of Lilly (1992) gives an equation for the model coefficient in the dynamic mixed model, 


where 


< Hij) > 

< 2 MuMki > 


L UiUj UjUj , 


Mjj = a 2 


D 


DijC’B — 


D 


Tlij C B ? 


Hij = ( UiUj - UiUj) - ( UiUj - UiUj), 

a is the ratio of the test and grid filters 

a = A/A, 


(2.46) 

(2.47) 

(2.48) 

(2.49) 

(2.50) 


and < > indicates averaging on horizontal planes. Equivalently, for the dynamic Smagorin- 
sky model, 


CA 2 = -- 


< MijCq > 


(2.51) 


< 2 MkiMki > 

The models outlined above use the original dynamic procedure of Germano et al. 
(1991) and require averaging in homogeneous directions to remain stable. Ghosal et al. 
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(1995) recast the dynamic procedure as a constrained variational problem and used this 
approach to develop a “dynamic localization model” that does not require averaging. 
This model, however, requires the solution of an integral equation and is computationally 
expensive. A less expensive alternative proposed by Piomelli & Liu (1995) involves finding 
an approximate solution to the integral equation by using the value of C at the previous 
time step to give a first approximation C* . Applying their procedure gives localised 
versions of the models. The coefficient for the localised dynamic Smagorinsky model is 
given by 


C = — 


(L 


*a 

ij 


- 2C*B. lj )A lJ 

2 AkiAki 


(2.52) 


where, 


Aij = a 2 A 2 


D 


It i \ ('n , 


Bij = A 2 


D 


DijCs- 


The coefficient for the localised dynamic mixed model is 

(£« - Hij - 2 C%)A y - 


C = — - 


O A. . A.. 


(2.53) 

(2.54) 


(2.55) 


with Aij and B t j given by Eqs.(2.53) and (2.54) respectively and H t j given by Eq.(2.49). 

We note here that by using a horizontal test filter we have reduced the complexity of 
our models considerably. With three-dimensional test filtering the models would be of 
similar complexity to the equivalent dynamic models for compressible flow whereas, when 
a horizontal filter is used, the models for the anelastic equations become similar to those 
for incompressible flow. In fact, if the discretisation operator is taken as being equivalent 
to the first level of grid scale Favre filtering (that is, we assume implicit filtering at the 
grid scale), then the only difference in the implementation occurs in the mixed model, 
where the second grid scale filter in the formula for Hjj (Eq.(2.49)) is a Favre filter that 
must be applied explicitly. 

An interesting feature of the localised dynamic model is that, if unconstrained, it can 
calculate negative values for the model coefficient, and thus predict regions of negative 
eddy viscosity. While such regions may be interpreted as regions of backscatter, or energy 
transfer from small to large scales, Ghosal et al. (1995) and Carati et al. (1995) note that 
the dynamic Smagorinsky model has no information about how much energy is actually 
available in the subgrid scales. Consequently, there is no mechanism by which such a 
“counter-gradient diffusion” model for backscatter can be turned off once all the subgrid 
scale energy has been removed. For this reason we chose to clip the eddy viscosity at 
zero. For the localised model backscatter is therefore included. For the plane averaged 
dynamic Smagorinsky model, however, the effect of backscatter is included in an average 
sense through a reduction of the plane averaged model coefficient. As noted above, in the 
mixed model the Leonard term allows backscatter, so there is no need to resort to using 
a negative eddy viscosity. 

The Smagorinsky and mixed models for the subgrid scale fluxes of the scalar variables 
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are written 


7 <t> — ~ £oC^A 2 
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o ^ B •> 
OX a 


7 <t> = Qo(<t> u j - <jwj) - QqC^/S. 2 
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3</> r 

W t/B- 
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(2.56) 

(2.57) 


Here <f> is a generic scalar representing or q t . The mixed model includes a term repre- 
senting the scalar equivalent of the modified Leonard term, 


L %j = eo{(t>Uj -<t>Uj). (2.58) 

Equations for the model coefficients are then written in a form analogous to that used 
for the subgrid scale stress. For the plane averaged dynamic Smagorinsky model, 


^ . o < FjE-j > 

0 “ ~<F k F k >’ 


where the resolved flux Ej is given by 
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For the dynamic mixed model, 
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where 


Gj = (< inij — <t>Uj) — (<t>Uj — (/>Uj). 

The coefficient for the localised dynamic Smagorinsky model is given by 
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The coefficient for the localised dynamic mixed model is 

{Ej - Gj - ClQj)Pj 


CU = — 


PkPk 


(2.59) 

(2.60) 
(2.61) 
(2.62) 

(2.63) 

(2.64) 

(2.65) 

(2.66) 

(2.67) 


with Pj and Qj given by Eqs.(2.65) and (2.66) respectively, and Gj given by Eq.(2.63). 

Close to the rough surface at the bottom of the domain the size of the turbulent eddies 
decreases and the turbulence models described above do not give the correct dissipation. 
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The standard Smagorinsky model is over-dissipative while the dynamic models are under- 
dissipative. With the standard Smagorinsky model we use a damping function suggested 
by Mason (1994), 

A^2 = A2 + ( k(z + z 0 )) 2 ’ (2 ' 68) 

where k is the Von Karman constant and z o the roughness length. With all models we 
also add a surface layer stress term suggested by Brown et al. (2001). This stress takes 
the form 


Tsfc 


C s f c a(z)p 0 \u\4>dz, 


(2.69) 


where C s /c is a scaling factor, a(z) a smoothing function, it the horizontal wind speed 
and </> the dependent variable. At present we follow Chow & Street (2002) and choose 
C s f c to give T s f c equal to half the surface flux at the first grid point above the wall. 


3. Simulation results 

The test case chosen for the simulations is the first research flight (RF01) of the second 
Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) field experiment. The 
computational domain extends 3.2km in the two horizontal directions and 1.5km in the 
vertical direction, with 96 x 96 x 128 cells in the x, y and 2 directions respectively. Cell- 
spacing is uniform in the horizontal directions and stretched in the vertical direction to 
give cells of height 5m close to the bottom surface and in the vicinity of the inversion. 
This grid is typical of that currently used in large eddy simulations of the atmospheric 
boundary layer. The spatial discretisation uses a second-order centered scheme with 
no flux limiting, while time integration uses a second order Runge-Kutta scheme. The 
simulation is integrated over a period of four hours. Further details of the simulation 
set-up are described in Kirkpatrick et al. (2003) and will not be repeated here. 

In the following we compare results obtained with four of the SGS turbulence models 
described above: the standard Smagorinsky (SM), dynamic Smagorinsky (DSM), dy- 
namic mixed (DMM) and localised dynamic Smagorinsky model (LDSM). The standard 
Smagorinsky model uses a constant coefficient of Cs = C x !' 2 = 0.18 and a Prandtl 
number of Pr sgs = 0.4. For reference, we also show results obtained with no SGS model. 

Figure 1 shows vertical profiles of condensed water, liquid water potential temperature, 
eddy viscosity and diffusivities and SGS Prandtl number predicted by the various models. 
These profiles are calculated by averaging over horizontal planes as well as over a time 
period of 30 minutes. All profiles are calculated for the last half hour of the simulation, 
that is for t = 3.5 — 4h. The dynamic Smagorinsky models give quantitatively similar 
profiles to the standard Smagorinsky model. LDSM gives somewhat higher average eddy 
viscosity and diffusivities than DSM. This is because negative values of eddy viscosity 
predicted by LDSM are clipped, whereas with DSM conditions that would produce lo- 
cal negative values have the effect of reducing the plane average. The unclipped model 
coefficient profiles for DSM and LDSM (not shown) were found to be almost identical. 
The dynamic mixed model gives lower values of eddy viscosity and diffusivity because 
the Leonard term component of the SGS stress is included separately. All models predict 
a SGS turbulent Prandtl number close to the standard value of 0.4 through most the 
boundary layer, although DSM predicts an increase to a value of approximately 2 in the 
region just below the inversion. The reason for this behaviour, which is not seen with the 
other models, is not clear at present. 










LES of the cloud-topped atmospheric boundary layer 


33 


880 

870 


a 

860 

e 

850 


840 

1.05 

o 1.00 
*-+^ 

O 

d 

* 

"O 

.2 0.95 


0.90 
100 

_ 80 

1 * 

a 

jl 60 

0. 

5j 40 
20 
0 

0 12 3 4 

Time (h) 

Figure 2. Time series of inversion height, cloud fraction and liquid water path: no SGS 

model, ■ • • SM, DSM, LDSM, DMM. 

Figure 2 shows time series of inversion height Zi nv , liquid water path LWP and cloud 
fraction. The entrainment rate E can be estimated from the rate of change of Zi nv using 
the formula, 

/7 2, ■ 

T7 1 n /o i \ 

E — dt ^sub’ (3.1) 

Here w su b is the subsidence velocity, which for the present case is estimated to be 
— 0.3cms _1 . dz inv /dt is calculated as the average rate of change over the period t = 2— 4h. 
The dynamic Smagorinsky models give an entrainment rate of approximately 0.36- 
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Figure 3. Eddy diffusivity, buoyancy frequency, and SGS and resolved buoyancy production 
at cloud-top: • • • SM, DSM, LDSM, DMM. 


0.37cms _1 , while the dynamic mixed model gives 0.40cms _1 . These results are in excellent 
agreement with the entrainment rate E = 0.38±0.04cms _1 estimated from observational 
measurements. The standard Smagorinsky model, however, overpredicts entrainment sig- 
nificantly giving E w 0.55cms _1 . The entrainment rate with no SGS model is 0.34cms _1 , 
which is slightly lower than that obtained with the dynamic models, but still within the 
error bounds of the observational measurements. The observed liquid water path and 
cloud fraction remained approximately constant during the experiment. This is also well 
predicted with dynamic models and no SGS model, while with the standard Smagorinsky 
model the results for these parameters again deviate significantly. 

To help explain these differences, Figure 3 shows a comparison of the eddy diffusivity 
Khi squared buoyancy frequency N 2 , resolved buoyancy production B and subgrid scale 
buoyancy production B sgs . Subgrid scale buoyancy production is approximated as 

B sgs « -K h N 2 . (3.2) 

The buoyancy frequency is calculated using a method described by Stevens et al. (2002). 
This approach is an extension of the work of MacVean & Mason (1990) and takes into 
account changes in stability that occur when subgrid scale mixing causes condensation or 
evaporation. Figure 3 shows the eddy diffusivity predicted by the standard Smagorinsky 
model decreasing through the inversion as a result of the stable stratification which drives 
the buoyancy correction coefficient Cb in the model to zero (see Eq.(2.16)). There is, 
however, a significant overlap region in which both the eddy diffusivity and buoyancy 
frequency are large. From Figure 3 and Eq.(2.16) it is clear that this overlap region 
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Figure 4. Second and third moments of vertical velocity. Simulation results (resolved component 

only): no SGS model, • • • SM, DSM, — • — LDSM, — • • • — DMM. Measurements: 

radar data (open circles), in situ measurements (filled circles). 


also constitutes a zone of negative subgrid scale buoyancy production. Meanwhile the 
resolved buoyancy production with SM is similar to that seen with the other models. 
Negative buoyancy production is associated with conversion of mechanical energy into 
potential energy. In the present case this is interpreted as turbulence doing work to 
mix the warmer air above the inversion with cooler boundary-layer air below. Thus, the 
extra —B sgs provided by the standard Smagorinsky model is the cause of the enhanced 
entrainment rate seen with this model. Excessive entrainment of warm, dry air across the 
inversion then leads to the observed reduction of cloud fraction and liquid water path. 

The same effect is not seen with the dynamic models which damp the model coefficient 
so that it becomes zero just below the inversion. Consequently there is no overlap region 
and no region of significant negative SGS buoyancy production. By coincidence, the 
same effect is also achieved when no SGS model is used, which explains the similarity 
between the results obtained with the dynamic models and no SGS model for large scale 
parameters such as entrainment rate, liquid water path and cloud fraction. 

Figure 4 shows vertical profiles for the second and third moments of the resolved ver- 
tical velocity component compared with the observational measurements. The second 
moment is well predicted in the dynamic model simulations, with numerical results in 
excellent agreement with observations. The shape of the profiles obtained with the dy- 
namic models is also typical of a well-mixed boundary layer as is expected for nocturnal 
stratocumulus. In contrast, the simulation using the standard Smagorinsky model does 
not agree as well with the observations and gives a two-peaked profile. Such a profile 
indicates a decoupling of the cloud layer and sub-cloud layer which is not apparent in 
the observational results. 
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The differences in profiles of the third moment are even more striking. The standard 
Smagorinsky model gives reasonable results in the surface layer, but does not give the 
negative values observed in the cloud layer. The dynamic models, on the other hand, give 
good agreement with the observations throughout the boundary layer, although, close to 
the surface they appear to under-predict the magnitude of this statistic somewhat. More 
work is required on the surface layer model to improve this. 


4. Future work 

The results presented in this paper clearly demonstrate the utility of dynamic turbu- 
lence models for atmospheric simulations. There is, however, much work still to be done. 
At present we are using a second-order accurate scheme for the spatial-discretisation. 
Consequently the numerical errors are of a similar magnitude to the terms associated 
with the subgrid scale models. The authors are currently working on a fourth-order 
scheme to resolve this issue. More work is also required in refining the surface layer 
model, especially where it is applied in combination with the dynamic models. As dis- 
cussed above, dynamic models tend to underpredict the dissipation in the surface layer. 
This is because as the surface is approached the energetic eddies are no longer resolved 
and the assumptions underlying the dynamic procedure break down. While the current 
surface layer model helps to overcome these issues, the surface layer velocity and scalar 
profiles still deviate significantly from profiles predicted by similarity theory. 
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